load IntnlData MV;
MV(MV<=0) = NaN;
NumPort = 10; MinStocks = 50;

TFirst = find(ym==199303);
D = [];
for cntry = [55 50 51]
    RPort = NaN(NumPort+1,2,T);
    NPort = NaN(NumPort+1,T);
    for t = TFirst-2:T-2
        R11_t = Return11(:,t);
        MV_t1 = MV(:,t+1);
        idx0 = isfinite(R11_t) & isfinite(MV_t1) & CountryRegionDummy(:,cntry) & ~IsREIT & ~IsMicro(:,t,cntry);
        if sum(idx0) < MinStocks*NumPort %number of stocks at this stage should be enough to populate all portfolios
            continue;
        end

        R11_t(idx0) = winsorize(R11_t(idx0),WinsorLevelX);
        RNext = Return(:,t+2); 
        if WinsorizeYRet
            RNext(idx0) = winsorize(RNext(idx0),WinsorLevelY);
        end
        MVNext = MV(:,t+1);    MVNext(idx0) = winsorize(MVNext(idx0),WinsorLevelY);
        
        % neutralize country effects
        for cc = 1:49
            idx_sort_cc = CountryRegionDummy(:,cc) & idx0;
            if sum(idx_sort_cc) > 0
                R11_t(idx_sort_cc) = R11_t(idx_sort_cc) - mean(R11_t(idx_sort_cc),'omitnan');
            end
        end
        
        Rprev_sort = sort(R11_t(idx0));
        NinPort = length(Rprev_sort)/NumPort;
        BrkPts_n = Rprev_sort([1 floor(NinPort*(1:NumPort))]); BrkPts_n(end) = Inf;
        
        for n = 1:NumPort
            idx_n = R11_t>=BrkPts_n(n) & R11_t<BrkPts_n(n+1) & idx0;
            if sum(idx_n) >= MinStocks %ensure that number of stocks in this portfolio is enough
                ret = RNext(idx_n); mv = MVNext(idx_n);
                idx = isfinite(ret) & isfinite(mv);
                ret = ret(idx); mv = mv(idx);
                if sum(idx) >= MinStocks % ensure again that we have returns for enough stocks
                    RPort(n,1,t+2) = sum(ret.*mv)/sum(mv);
                    RPort(n,2,t+2) = mean(ret);
                    NPort(n,t+2) = sum(idx);
                end
            end
        end
    end
    % take data only where corner portfolios are populated
    tmp = [RPort(1,1,:) RPort(NumPort,1,:)];
    idx = all(isfinite(tmp),2);
    RPort(:,:,~idx) = NaN; NPort(:,~idx) = NaN;
    RPort(NumPort+1,:,:) = RPort(NumPort,:,:) - RPort(1,:,:);
    
    Y = squeeze(RPort(NumPort+1,1,:)); Y(~isfinite(Y)) = [];
    D1 = [100*mean(Y); 100*median(Y); 100*std(Y); skewness(Y); kurtosis(Y); 100*min(Y); 100*max(Y)];
    D = [D D1];
end
writematrix(D,FileName, 'Sheet', SheetName, 'Range','b2:d8', 'AutoFitWidth',false);
